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Abstract 

Background: The circadian clock is a critical regulator of biological functions controlling behavioral, physiological 
and biochemical processes. Because the liver is the primary regulator of metabolites within the mammalian body 
and the disruption of circadian rhythms in liver is associated with severe illness, circadian regulators would play a 
strong role in maintaining liver function. However, the regulatory structure that governs circadian dynamics within 
the liver at a transcriptional level remains unknown. To explore this aspect, we analyzed hepatic transcriptional 
dynamics in Sprague-Dawley rats over a period of 24 hours to assess the genome-wide responses. 

Results: Using an unsupervised consensus clustering method, we identified four major gene expression clusters, 
corresponding to central carbon and nitrogen metabolism, membrane integrity, immune function, and DNA repair, 
all of which have dynamics which suggest regulation in a circadian manner. With the assumption that transcription 
factors (TFs) that are differentially expressed and contain CLOCK:BMAL1 binding sites on their proximal promoters 
are likely to be clock-controlled TFs, we were able to use promoter analysis to putatively identify additional clock-controlled 
TFs besides PARF and RORA families. These TFs are both functionally and temporally related to the clusters they regulate. 
Furthermore, we also identified significant sets of clock TFs that are potentially transcriptional regulators of gene clusters. 

Conclusions: All together, we were able to propose a regulatory structure for circadian regulation which represents 
alternative paths for circadian control of different functions within the liver. Our prediction has been affirmed by functional 
and temporal analyses which are able to extend for similar studies. 

Keywords: Circadian rhythm, Microarray analysis, Gene expression, Consensus clustering, Promoter analysis, Transcription 
factor, Circadian regulation 



Background 

The circadian clock is one of the most critical biological 
regulators for all living entities controlling behavioral, 
physiological and biochemical processes [1]. Acting like 
a multifunctional timer with a roughly 24-h cycle [2], 
the circadian clock offers fitness advantages by endowing 
organisms with anticipatory mechanisms predicting peri- 
odic events, as opposed to responding in a reactive way to 
external signals [3]. Robust circadian regulation is associ- 
ated with fitness, while abnormal rhythms which are char- 
acteristic of stress are often linked with illness [4,5]. As an 
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example, metastatic colorectal cancer patients with severe 
alterations of the rest-activity circadian rhythm were pre- 
dicted to have a fivefold increase in the risk of death com- 
pared to a normal rest- activity pattern [6]. 

In mammals, the circadian clock is composed of auto- 
regulated transcription-translation feedback loops whose 
primary elements are transcription factors (TFs) includ- 
ing the basic helix-loop-helix transcription factor CLOCK 
and the aryl hydrocarbon receptor nuclear translocator- 
like BMAL1, which drive the expression of period (PER) 
and cryptochrome (CRY). PER and CRY form the negative 
limb of the feedback loop, repressing the expression of 
CLOCK-and BMALl-induced transcription [7]. Evidence 
suggests that circadian gene expression is achieved 
through a transcriptional cascade in which the expres- 
sion of other TFs responsible for the majority of rhythmic 
tissue-specific expression are primarily driven directly by 
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the main clock transcription factor complex CLOCK: 
BMAL1 [8]. Recently, this hypothesis has been reinforced 
by experiments showing that TFs whose expression pat- 
terns exhibit daily oscillation are direct targets of the 
CLOCK:BMALl complex [9,10]. However, those studies 
only explored the transcriptome with one main clock 
transcription factor mutated and thus the entire picture of 
how the system behaves or how target genes are con- 
trolled is still a question. Consequently, among the aims 
of our study is to gain insights into the broader transcrip- 
tional regulation of circadian rhythms. 

Liver is one of the most important organs regulating 
the metabolic activity of the body as well as producing 
acute phase proteins (APP, class of proteins responding 
to inflammatory stressors) in response to external stress 
[11]. Genome-wide transcriptome-mapping studies have 
revealed that the fraction of diurnally regulated liver tran- 
scripts could amount up to 10%, with genes exhibiting cir- 
cadian rhythms participate in various functions e.g. cell 
signaling, energy metabolism, amino acid, lipid and choles- 
terol metabolism, carbohydrate transport and metabolism, 
DNA replication, protein synthesis, signal transduction 
mechanisms, and immune response [8,12-15]. Due to the 
significance of the circadian periodicity and the severe out- 
come of disorders following its disruption, a global per- 
spective of the transcriptional dynamics and regulatory 
structure of circadian genes in liver may offer a molecular 
framework for studies on the hepatic circadian rhythm 
gene regulation. 

In this study, we address two critical questions including 

(1) whether we can identify and functionally characterize 
circadian-regulated genes in an unsupervised manner, and 

(2) if the potential transcriptional regulatory mechanisms 
of circadian rhythms can be hypothesized. In order to ob- 
tain the necessary samples for analysis, animals were sacri- 
ficed at 9 am, 11 am, 1 pm, 5 pm, 1 am, and 9 am (~0, 2, 
4, 8, 16, and 24 hr zeitgeber time) via exsanguinations and 
livers were harvested, allowing for the acquisition of tran- 
scriptional data from microarray readouts for those time 
points. Utilizing our previous consensus clustering [16] 
and techniques in promoter analysis [17], our results actu- 
ally demonstrate that clusters of co-expressed genes do 
vary in a circadian manner and contain functionalities 
consistent with previous studies including immune func- 
tion, cell repair, metabolism, and DNA repair. Additional 
clock-controlled TFs are also predicted. Interestingly, we 
observed that many clock- controlled TFs relevant to the 
regulation of genes within a circadian pattern show similar 
transcriptional dynamics. Furthermore, the predicted 
regulatory network and regulatory motifs of circadian 
patterns are shown to be consistent with the biological 
significance of the corresponding gene set and the cir- 
cadian time of the transcriptional pattern. Finally, the 
identified transcription factors are intimately related to 



innate immunity and response to infection, providing 
further evidence that circadian rhythms play an im- 
portant role in recovery from injury. While it is under- 
stood that recent evidence points to the important role 
alternative processes play for driving the induction of 
oscillations [18-20], the proposed approach presents a 
generalizable framework for assessing the regulatory struc- 
ture underlying the transcriptional sources of oscillations. 

Results 

Circadian dynamics and functions 

Male Sprague-Dawley rats (Charles River Labs, Wilmington, 
MA) were housed in a temperature-controlled environment 
(25°C) with a 12-hour light-dark cycle and provided water 
and standard chow ad libitum. Because the animals are noc- 
turnal, this establishes a 12 hour LD cycle for the animals, 
where ZT0-ZT12 represents the resting portion of the day, 
while ZT12-ZT0 represents the active portion of the day. 
These animals were sacrificed at 9 am, 11 am, 1 pm, 5 pm, 
1 am, and 9 am (~0, 2, 4, 8, 16, and 24 hr zeitgeber time) 
and liver tissues were collected and frozen for microarray 
analysis (n = 3 per time point per group). The tissues were 
lysed and homogenized using Trizol, and the RNAs were 
further purified and treated with DNase using RNeasy col- 
umns (Qiagen). Then cRNAs prepared from the RNAs of 
liver tissues using protocols provided by Affymetrix were uti- 
lized to hybridize Rat Genome 230 2.0 Array (GeneChip, 
Affymetrix) comprised of 31,099 probe sets. In order to 
analyze the microarray data, we first applied the ANOVA 
test to filter for differentially expressed genes (p-value< 
0.001, q-value < 0.01). 545 probe sets are identified as differ- 
entially expressed. Given the assumption that a significant 
set of co-expressed genes can characterize some crit- 
ical transcriptional dynamics and biological signifi- 
cance, we explored the potential of our prior work to 
identify highly co-expressed gene clusters [16]. Specif- 
ically, the algorithm performs a consensus clustering 
and a trivial cluster removal procedure resulting in 
four significant clusters composed of 153, 64, 52, and 
83 probesets (or 132, 42, 46, 70 genes respectively- 
some probesets are unmapped) that characterize four 
critical circadian dynamics of the homeostatic system 
in rat liver (Figure 1). The gene list and expression 
levels for each probe set are provided in the Additional 
file 1. Additionally, in order to test whether our identi- 
fied genes express in a circadian manner, we applied 
the Fisher s exact g-test which is discussed extensively 
in Wicher et al. [21] to assess the significance of the 
periodic transcription patterns. From the original ex- 
pression vector of each probeset with 3 replicates at 
each time-point, a new vector is created for the test 
where replicates at time-point t (hr) are randomly dis- 
tributed to three consecutive time-points t, t + 24, and 
t + 48. The result is that with q-value <0.01, all 352 
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(abat, agxt2, aldh4a1, gls2) 

- Glycine, serine and threonine metabolism 
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- beta-Alanine metabolism (abat, dpys) 
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- Energy production (ppargda, coqIOb) 

- Cell-cell junction (tand, epb41l5) 

- Membrane structure (ubr4, tmem106b) 



- Chemokine signaling (chuk, gnai3, stat3) 

- RIG-l-like receptor signaling (chuk, ddx3x) 

- Toll-like receptor signaling (chuk, myd88) 



- Nucleotide excision repair 
(Ioc1003621,rfc2,rfc5) 

- Mismatch repair (rfc2, rfc5) 

- DNA replication (rfc2, rfc5) 



Transcription Patterns 



0 4 8 12 16 20 24 
0 4 8 12 16 20 24 
0 4 8 12 16 20 24 




0 4 8 12 16 20 24 



Color key: 
Row z-score 



Time (hr) 



-1.5 0 1.5 

Figure 1 Circadian patterns in homeostatic rat liver. Left panel displays the heat map of expression of 153, 64, 52 and 83 probesets in 4 
patterns at 0, 2, 4, 8, 16, and 24 hr. Red color indicates the highest level while green indicates the lowest level of expression. Middle panel 
consists of enriched functions characterizing biological significance of each circadian pattern. Right panel shows the average normalized (z-score) 
expression profiles of probesets in four corresponding patterns. Error bars are two standard deviations of all probeset transcript levels at each 
time-point in each corresponding pattern. 



selected probesets in 4 clusters pass the test, implying 
that these 290 genes do reflect periodic expression pat- 
terns. Thus they are so-called as circadian- relevant genes 
in this context. 

The detailed dynamics of the four transcription pat- 
terns are shown in Figure 1 with the average expression 
profiles of all probesets in each corresponding cluster. 
Functional annotation was performed to develop a suc- 
cinct description of the key roles of each pattern. In 
brief, pattern 1 (132 genes) is composed of genes whose 
activity peaks in the morning (ZTO ~ 9:00 am) and is 
enriched by functions associated with metabolism, in- 
cluding (i) amino acid metabolism such as nitrogen me- 
tabolism (car 14, carSa, cth, gls2), alanine, aspartate and 
glutamate metabolism (abat, agxt2, gls2, loc641316), gly- 
cine, serine and threonine metabolism (agxt2, cth, gldc), 
arginine and proline metabolism (agmat,gls2), beta- 
alanine metabolism (abat, dpys), and (ii) carbohydrate 
metabolism such as glycosaminoglycan degradation (gns, 
gusb), propanoate metabolism (abat,acacb), starch and 
sucrose metabolism (gusb, pygl). In addition, many other 



genes related to primary metabolism and are known to 
oscillate in a circadian manner e.g. bmal, gcjr, dhtkdl, 
cth, gaba, gldc, kynu, akrlcl8 [22], chka and agxt [23]. 

Pattern 2 contains 42 genes whose expression levels 
reach the maximum during the day time (ZT2 ~ 
11:00 am) and appears to approach baseline before the 
dark phase. This pattern consists of genes with diverse 
functions including energy production (ppargda, coqIOb), 
cell-cell junction (tancl, epb4.1 IS), and membrane struc- 
ture (ubr4, tmeml06b). The genes associated with cell 
junction and component of membrane play a critical role 
in the integrity of the endothelial barrier function. 
Additionally, this pattern also contains several tran- 
scription factors involved in regulating circadian dy- 
namics e.g. CREM (a transcription factor involved in 
output clock function and melatonin synthesis in the 
mouse [24]), KLF9 [25]. In a similar manner, 46 genes 
in pattern 3 display a maximum activity during the 
light phase and are primarily composed of genes asso- 
ciated with the immune system. Genes in this temporal 
class are involved in the Chemokine signaling pathway 
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(chuk, gnai3, stat3), RIG-I-like receptor signaling path- 
way (chuk, ddx3x), and Toll-like receptor signaling 
pathway {chuk, myd88). Also, there are numerous genes 
that were shown to have circadian dynamics and related 
to immune function in previous studies such as btgl [26], 
sfpq [27], and f<bpS [22]. Finally, pattern 4 (70 genes) ex- 
hibits a reduced activity during the day time and reaches 
its peak activity around midnight. It is characterized by 
functions associated with DNA replication and repair, 
such as mismatch repair {rfc2, rfc5), DNA replication 
(rfc2, rfcS), and non-homologous end joining (rad50). Fur- 
thermore, timeless [28] and calml [26] are also known as 
genes with circadian rhythmicity that vary within this 
pattern. 

Putative transcriptional regulators controlling circadian 
dynamics 

In order to gain insights into the underlying regulatory 
mechanisms of circadian rhythms, we first analyzed the 
promoter regions of these putatively circadian-regulated 
genes in an attempt to predict putative transcriptional 
regulators and potential regulatory relationships. From 
290 genes (of 352 probesets) in four transcription pat- 
terns, we further analyzed 134 genes that have sufficient 
information of orthologous promoters (at least 3 ortho- 
logous promoters) (see Methods). Table 1 shows identi- 
fied putative transcriptional regulators associated with 
each circadian pattern. Since transcription factors are 
characterized by pleiotropic effects, it is reasonable to 
observe an overlap across the regulation of various tran- 
scription patterns [29]. While comparing these regula- 
tory combinations, we observe that some TF families 
seem to be common regulators for all rhythms, suggest- 
ing that circadian dynamics of these genes may be regu- 
lated by a common underlying regulatory mechanism. 
These TFs are shown in bold in Table 1. Since they are 
commonly present as transcriptional regulators across 
all four patterns, they may be involved in the control of 



dynamics of all circadian-relevant genes identified in the 
homeostatic system. On the other hand, clock TFs are 
hypothesized to be transcriptional regulators relevant to 
the regulation of circadian dynamics, however, because 
these common TFs may or may not be clock TFs, it is 
possible that genes encoding them are regulated by clock 
TFs and in turn they regulate their target genes to ex- 
press in a circadian manner. This proposed assumption 
is further tested in following sections. 

Transcriptional patterns of clock-controlled TFs and their 
temporal organization 

To further investigate the underlying regulatory mecha- 
nisms of circadian rhythms, we explored the hypothesis 
that CLOCK:BMALl, a core clock transcription factor 
family including CLOCK, BMAL1, NPAS2 [30], is respon- 
sible for the regulation of circadian rhythms through cis- 
regulatory modules by coordinating with other TFs whose 
dynamics are clock controlled. By assuming that TFs 
which are target genes of CLOCK:BMALl and also differ- 
entially expressed in homeostatic rat liver are potentially 
clock-controlled TFs (ccTFs), we scanned the list of 545 
differentially expressed probesets to identify additional 
clock-controlled TFs (see Methods). A list of 10 additional 
ccTFs (distributed in seven TF families including STAT, 
FKHD, CREB, RXRF, KLFS, AHRR and EBOX) was identi- 
fied. These transcription factors have dynamics that, 
when plotted in circadian zeitgeber time (ZT) where 
ZT0 (-9:00 am) is defined as the onset of diurnal ac- 
tivity, show an interesting progression of their peak ac- 
tivities. The plot of clock-controlled TFs vs. their peak 
circadian time and the individual profiles of the ccTFs 
are shown in Figure 2. Worth noting is the fact that 
many of these ccTFs exhibit similar transcriptional 
patterns. Although there appears to have no relation- 
ship between transcriptional levels and TF activity, this 
similarity may propose some cooperation in the tran- 
scriptional regulation of circadian-relevant genes. 



Table 1 Putative transcriptional regulators controlling circadian-relevant genes 

Circadian Putative transcriptional regulators 

patterns 

1 63 genes / 132 genes* 

AP1R, AP4R, CAAT, CLOX, CREB, E2FF, ETSF, EVI1, FKHD, GATA, GREF, HOXF, MYBL, MYT1, NKXH, NR2F, OCT1, PLAG, RUSH, RXRF, 
SORY, STAT, XBBF 

2 15 genes / 42 genes 

AHRR, CLOX, CREB, CTCF, E2FF, EREF, ETSF, EVI1, FKHD, GCMF, GLIF, HAND, HEAT, HOMF, HOXF, MYBL, NKXH, NR2F, OCT1, PARF, 
PAX5, PERO, RXRF, SORY, SP1F, XBBF 

3 24 genes / 46 genes 

AP1R, CREB, E2FF, EGRF, ETSF, GATA, HEAT, HOXF, IRFF, MYBL, NKXH, NR2F, PERO, RXRF, STAT, TBPF 

4 32 genes / 70 genes 

CREB, E2FF, ETSF, EVI1, NR2F, RXRF, SORY, SP1F 

*:number of genes with orthologous promoter information vs. the total number of genes in the group; TFs shown in bold are those present as regulators in at 
least three over four groups of genes. 



Nguyen et al. BMC Bioinformatics 2014, 15:83 
http://www.biomedcentral.eom/1 471 -21 05/1 5/83 



Page 5 of 14 



(a) 



T, STAT5A 



A FOX03 



TF Family 
RXRF 



(b) 



NR1I3 



0 4 8 12 16 20 24 



0 4 8 12 16 20 24 





0 4 8 12 16 20 24 



0 4 8 12 16 20 24 



KLF9 jl RARB 
-2 H 1 1 ^^^^^^^^^^ -2 H \ 1 ^^^^^^^^^^ 



CTAT 
O 1 Al 


STAT3 


CREB 




ATF6 




EBOX 






AHRR 


AHR 


RXRF 


RARB 


KLFS 


KLF9 


CREB 


CREM 




F0XN3 
F0X03 


FKHD 




STAT 


STAT5A 












► 



ZTO ZT2 



0 4 8 12 16 20 24 



t r~ 

0 4 8 12 16 20 24 



ZT4 ZT8 
Zeitgeber Time 



ZT16 




-1 - 



-1 - 



12 16 20 24 




STAT3 



12 16 20 24 



0 4 8 12 16 20 24 



J NR1I3 

0 4 8 12 16 20 24 

Figure 2 Transcription patterns of additional clock-controlled TFs. (a) Expression profiles of each TF in the homeostatic rat liver; the 
horizontal axis is the zeitgeber time (hr) and the vertical axis is the z-score normalized expression levels. Error bars are two standard deviations of 
expression values at each particular time-point, (b) The temporal organization of clock-controlled TFs. TFs with similar binding sites are combined 
into a TF family in the computational context. They are arranged following their corresponding peak circadian times. 



Putative clock TFs regulating the transcription of 
circadian-relevant genes 

As it has been hypothesized that TFs usually regulate 
gene expression in a combinatorial manner rather than 
in isolation [31,32], clock TFs, including the core clock 
TFs and clock-controlled TFs, are also hypothesized to 
cooperate to regulate the transcriptional levels of circadian- 
relevant genes. We thus searched for significant groups of 
clock TFs that are significantly over-represented on the cor- 
responding promoter set of genes in each pattern. From 
290 selected genes in four circadian patterns, 234 genes 
have promoter annotation in Genomatix database (see 
Supplementary for the gene list). 620 Rattus norvegicus's 
promoter sequences were extracted. For each circadian pat- 
tern, a list of common sets of clock TFs were identified and 
only those that are statistical significant compared to the 
GC-matched background are selected (see Methods) 
(Table 2). A union of these TF sets would provide a putative 
set of clock TFs associated with the regulation of corre- 
sponding genes characteristic of each circadian pattern. 



Table 2 Putative clock TFs involved in the transcriptional 
regulation of circadian-relevant genes 



Circadian 


Sets of regulators 


p-value | 


Sets of 


p-value | 


patterns 




q-value 


regulators 


q-value 


1 109/132 


CREB, RORA, 


0.0078 | 


FKHD, PARF, 


0.0377 | 


genes* 


STAT 


0.0616 


RXRF 


0.0877 




CREB, RORA, 


0.0116 | 


RORA, RXRF, 


0.0409 | 




RXRF 


0.0616 


STAT 


0.0877 




PARF, RORA 


0.0191 | 


CREB, FKHD, 


0.0473 | 






0.0753 


RORA 


0.0877 


2 27/42 


AHRR, CREB, 


0.0077 | 


EBOX, CLOCK: 


0.0303 | 


genes 


STAT 


0.0616 


BMAL1 


0.0877 




AHRR, CREB, EBOX 


0.0079 | 


AHRR, RXRF, 


0.0327 | 






0.0616 


STAT 


0.0877 




AHRR, CREB, KLFS, 


0.0130 | 


AHRR, KLFS, 


0.0482 | 




RXRF 


0.0616 


STAT 


0.0877 


3 39/46 


CREB, EBOX, 


0.0428 | 






genes 


RXRF, STAT 


0.0877 







4 59/70 
genes 



None 



*:number of genes with promoter information vs. the total number of genes in 
the group. 
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Additionally, we tested whether the regulation of the com- 
mon set of transcriptional regulators identified using the 
foot-printing techniques in Table 1 is influenced by clock 
TFs (indicated in boldin Table 1). Results show that two 
sets of clock TFs including KLFS_STAT (p-value < 0.037) 
and FKHD_KLFS (p-value < 0.041) are likely to be relevant 
to the regulation of the set of common regulators L. This 
implies that genes in circadian pattern 4 may be indirectly 
regulated by clock TFs. 

Discussion 

Identified circadian dynamics and biological significance 

The liver expresses a diverse set of genes and previous 
research indicated that many hepatic genes are under 
direct or indirect circadian control [22]. Despite the fact 
that unsupervised pattern identification methods were 
used to identify dominant co-expressed gene sets, our 
results show that selected genes oscillate in a circadian man- 
ner, as shown in Figure 1. Functional characterization 
reveals that enriched functions are associated with me- 
tabolism, energy production, immune system and DNA 
replication and repair. 

The liver plays an important role in the orchestration 
of metabolic functions such as the regulation of meta- 
bolic fuels, protein synthesis, iron and vitamin storage, 
as well as nutrient intake [33]. As might be expected, an- 
imals subjected to a light-dark schedule display pro- 
nounced rhythms in glycogen content, with a peak 
occurring late in the night, following the main period of 
food intake. This is also in agreement with genes identi- 
fied in pattern 1, where enriched metabolic functions, 
especially glucose metabolism, reach their peak during 
the dark phase. Since the dark phase represents the ac- 
tive period for rats, where they are moving and eating, 
this increase in metabolism is indicative of both energy 
burn to fuel muscle usage, as well as digestive effects on 
the liver, and an influx of glucose and other amino acids 
to be processed through glycogen production and the 
urea cycle. Because the current study does not differenti- 
ate from fed animals and fasted animals, it is impossible 
to determine how much increased metabolic activity is 
due to demand (movement), versus supply (digestion). 
In addition, it has been suggested that the expression 
patterns of circadian rhythms in peripheral organs such 
as the liver are the result of delayed feedback loops of 
gene products from circadian regulators, suggesting that 
upstream regulation in addition to direct transcription 
factor binding can occur up to 6 hours prior to the peak 
of gene expression [34]. 

Energy production is also time-of-day dependent. It is 
reported that the production of ATP in rat liver is also 
regulated by the circadian rhythm and is lower in the 
dark phase [35]. Membrane integrity and cell to cell 
junctions are another critical circadian function that has 



been previously hypothesized to be important during 
rest [36], and contain significant circadian regulation 
[25]. Circadian pattern 2 is made up of a diverse array of 
genes whose functions are related to energy production 
and membrane structure, but are up regulated during 
the day. Many of these metabolic genes are associated 
with mitochondrial biogenesis specifically, which is a 
process that has been previously shown to occur during 
the resting phase of the circadian cycle [37]. As rats are 
nocturnal, their day cycle is indicative of their sleep 
cycle, this is in prior agreement and it suggests that 
mitochondrial biogenesis is a major liver function during 
the rest portion of their circadian rhythm. Genes that 
are associated with apoptosis are also found within this 
pattern such asbclafl, rarb, and dapkl and are further 
indication of cellular repair-cells which are too damaged 
to be simply repaired are killed and replaced in order to 
prepare the host for another day of activity. 

Immune system functions, including lymphocyte pro- 
liferation, natural killer cell activity, humoral immune re- 
sponse, absolute and relative numbers of circulating 
white blood cells and their subsets, cytokine levels, and 
serum Cortisol, are known to be entrained by circadian 
cues [27]. Also, it has been previously reported that 
TNFa secretion was significantly increased in burn vs. 
sham isolated splenocytes cells only when injury took 
place in the morning [38] further verifying that the re- 
sponse to injury is time-of-day dependent as a result of 
the circadian variability of the immune functions. This is 
confirmed by the selected gene set exhibiting pattern 3 
including&£gi, sfpq, fl<bpS, which is primarily related to 
immune function, and has consistent dynamics with pre- 
vious works [22,26,27]. However, further analysis of the 
dynamics observes a peak in immune function at the 
end of the day (ZT8 ~ 5:00 pm), which is the time rats 
are beginning to awaken. Over the course of the night, 
when the rats are active, immune function gradually de- 
creases, until it reaches a minimum when the night 
ends, and the rats rest. This implies that there may be a 
recovery mechanism associated with sleep that allows 
for the rejuvenation of immune function in order to 
ward off infection during active periods. 

The strategy of adaptive circadian clocks could be also 
timing of UV-sensitive cellular processes to occur at 
night to avoid UV-induced damage [39]. Our analysis 
showed that functions related to DNA replication/repair 
(pattern 4) achieve a maximum activity around mid- 
night, which is consistent with previous works [28]. 
However, unlike the previous patterns, the circadian dy- 
namics of these genes are not dependant on the habits 
of the host animal, but rather upon external conditions 
related to the environment. Unlike the other patterns, 
the genes contained in pattern 4 do not map to specific 
transcription factors directly, which suggests that this 
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group is under a different form of control than the 
others. The expression pattern of pattern 4 is similar to 
the late cryl phase of the circadian rhythm [40], which 
were identified to be under the control of newly identi- 
fied D boxes that overlap with E-boxes to provide delayed 
expression of genes which contributes to rhythmicity. Fur- 
thermore, the cryl expression pattern, whose phase re- 
sembles cluster 4, is heavily dependent on other circadian 
regulators for its robustness [34], which suggests that pat- 
tern 4 may be controlled through upstream transcription 
factor regulators through their gene products, in addition 
to direct TF binding. 

Putative regulatory motifs of circadian-relevant genes 

In order to investigate the underlying regulation behind 
the co-expressed circadian genes, promoter analysis was 
applied to identify key transcriptional regulators associ- 
ated with the regulation of each pattern. The initial as- 
sumption for this analysis is that core-clocks TFs are 
responsible for the regulation of circadian rhythms 
through coordination with other TFs whose dynamics 
are clock controlled (ccTFs). However, besides clock TFs 
(including core-clock TFs and ccTFs) there could be 
other TFs also involved in the regulation of transcrip- 
tional patterns of these circadian-relevant genes. As 
such, we first explored the underlying hypothesis in 
comparative genomics to predict relevant transcriptional 
regulators using the foot-printing technique (Table 1) 
[41,42]. Subsequently, we explored the initial assumption 
above by searching for statistical significant sets of clock 
TFs that are over-represented on the corresponding pro- 
moter set of genes in each pattern (see Methods). PARF 
and RORA are also included in this search because they 
are widely known to be regulated by core- clock TFs 
(Table 2-the 3 left rectangles in Figure 3a). Furthermore, 
to exhaustively examine if any set of clock TFs relevant 
to the regulation of common regulators of circadian 
genes, we applied the latter methodology to the set of 
common regulators identified in Table 1 (TFs in the cir- 
cle of Figure 3a). Results showed that a set of three 
clock-controlled TFs (STAT, FKHD, and KLFS) are asso- 
ciated with their regulation. These results are integrated 
and displayed in Figure 3a, which consists of TFs that 
are regulated by CLOCK:BMALl and directly regulate 
individual circadian patterns, as well as those that are 
regulated by ccTFs and affect the dynamics of all four 
patterns. 

Simplified networks are shown in Figure 3b, which sim- 
ply condense the transcription factors into three categor- 
ies: core clock transcription factors (CLOCK:BMALl), 
clock controlled transcription factors (ccTFs) and other 
transcription factors. Genes which exhibit transcription 
pattern 2 show the most circadian influence, with the cor- 
responding promoter set containing enriched binding sites 



for the core clock TFs, ccTFs, and other TFs. Meanwhile 
genes following pattern 4 has the least circadian influence, 
being only controlled indirectly through two levels of 
regulation by the core clock transcription factors, likely 
through their gene products [40]. Because of the overlap- 
ping levels of regulation associated with pattern 2, it can 
be expected that disruptions to the circadian cycle will be 
less effective at disrupting its gene expression, as those 
layers of regulation could potentially provide redundancy 
that preserves circadian gene expression. Meanwhile, pat- 
tern 4, can be putatively said to be more sensitive to circa- 
dian disruptions, due to the lack of backup regulation. 
This is supported by the fact that pattern 4 has a similar 
expression to the cryl gene, which is a core circadian 
regulator known for its sensitivity to changes in the overall 
circadian rhythm. Further, since pattern 2 is representative 
of host repair during rest, it makes sense that it would be 
protected against disruption: proper sleep is critical to 
health and disease recovery [43], and the redundancy of 
the host repair genes may reflect that need. Similarly, pat- 
tern 4 is associated with DNA repair, which has been 
shown to be severely affected following injury due to the 
increased presence of reactive oxygen species [44,45], 
known to damage DNA. Thus the putative sensitivity of 
this pattern to circadian disruption may be a result of al- 
tered DNA repair patterns during periods of stress. Genes 
in pattern 1 and 3 are both under similar circadian regula- 
tion, despite being significantly separated in circadian 
time, as discussed previously. These patterns primarily re- 
late to metabolism and innate immunity, two functions 
which have been shown to be closely linked through a 
phenomenon known as hypermetabolism [46]. These 
functions show significant circadian sensitivity, and be- 
cause of their peaks at different circadian times (ZT0 and 
ZT8), it is possible that due to different innate immune 
and metabolic states over the course of the circadian day, 
the time of injury can significantly impact the severity of 
the immune response and the onset of hypermetabolism. 
The hypothesis that circadian rhythms strongly impact 
immune responses following injury has been postulated 
before [47], but further experiments are needed in order 
to conclusively link the functions, and dynamics of genes 
in pattern 1 & 3 to injury responsive elements. 

Dynamic and functional agreement between clock- 
controlled TFs and their associated gene sets 

In Figure 2, it is worth noting that the transcription fac- 
tors that have early circadian time peaks (ZT0 to ZT4), 
namely STAT5A, FOX03, FOXN3, CREM, and RARB 
are all associated with the regulation of genes in tran- 
scription pattern 1, whose genes also peak very early in 
circadian time. Transcription factors which peak at 
somewhat later times (ZT4-ZT8) include RARB, KLF9, 
AHR, ATF6, and STAT3, which have all been associated 



Nguyen et al. BMC Bioinformatics 2014, 15:83 
http://www.biomedcentral.eom/1 471 -21 05/1 5/83 



Page 8 of 14 



CLOCK: BMAL1 



CREB, FKHD, 
PARF, RORA, 
RXRF, STAT 




0 4 8 12 16 20 24 0 4 8 12 16 20 24 



4 8 12 16 20 24 0 4 8 12 16 20 24 



(b) 




circadian- 
relevant genes 



CL0CK:BMAL1 



(1&3) 



ccTFs 




other 
TFs 



circadian- 
relevant genes 



CL0CK:BMAL1 
(4) 

ccTFs 

other 
TFs 



circadian- 
relevant genes 



Figure 3 Putative regulatory program of circadian-relevant genes (a) and proposed regulatory motifs (b) ccTFs are clock-controlled TFs. 

Regulation of rhythm 2 is corresponding with the first regulatory motif rhythm 1&3 with the second and rhythm4 with the last one. 



with pattern 2. This pattern peaks at a later time period 
than the first pattern, and this is reflected in the shift in 
circadian time peaks in its transcriptional regulators. Fi- 
nally, the transcription factors which peak at the latest 
circadian times (ZT8-Z16) include ATF6, STAT3, and 
NR1I3, which are all associated with pattern 3. This pat- 
tern peaks at the latest time point following the diurnal 



onset, and again, this is reflected in the shifted circadian 
time peaks of its regulators. In general, many of tran- 
scription factors control the same set of genes exhibiting 
a circadian pattern express a very similar circadian 
profile. 

Though many transcription factors have been identified, 
the three clock-controlled TF families which control the 
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functions of all of the patterns through binding sites on 
non-clock controlled TFs are FKHD, KLFS, and STAT 
(Figure 3). The FKHD transcription factor family includes 
the Fox-O transcription factor, which is well known for its 
protective role against oxidative stress within the body 
[48]. Similarly, KLFS is well known for its role in the regu- 
lation of anti polymicrobial bacteria activity within the 
body and metabolism in liver [25,49], while STAT has 
already been identified as a critical transcription factor 
during the outcome of sepsis [50]. In particular, KLF has 
been strongly associated with liver metabolism in particu- 
lar [25], which is reinforced by its regulation of cluster 2, 
which has many metabolic genes associated with it. Be- 
cause KLF is strongly associated with CLOCK and BMAL, 
its presence as a regulator of these oscillatory gene clus- 
ters suggests significant circadian regulation of these func- 
tions. Together these three mediators and their functions 
indicate that the overarching control of clock TFs to the 
patterns indirectly through non clock-controlled TFs is 
controlled primarily by the innate immune response, and 
may serve to coordinate the metabolic functions and toll 
like receptor functions found in the patterns in the event 
of an innate immune response. 

A TF family that is considered as a common regulator 
of genes across all four circadian patterns is CREB. It 
has been previously recognized to be prototypical stimu- 
lus induced in function [51], and therefore plays a wide 
variety of roles in different tissues. However, it has been 
strongly linked to functional circadian regulation through- 
out the body, and thus it is not surprising that CREB regu- 
latory elements are found in all four patterns, both 
directly and indirectly regulated by CLOCK:BMALl. 
Unique clock controlled TF families that regulate genes in 
pattern 1 are PARF and RORA. PARF includes well known 
clock TFs (e.g. DPB, TEF) which are directly targets of 
CLOCK:BMALl [52,53]. RORA is another well known 
circadian transcription factor family, which has also been 
shown to play a role in metabolism in many tissues, in- 
cluding the muscle and liver systems [54,55]. Additionally, 
RXRF transcription factors represent a family of nuclear 
receptors that target a variety of signaling pathways, which 
are primarily related to metabolism, cell differentiation, 
and cell death [56]. Thus, much like others associated with 
pattern 1, RXRF factor family has significantly overlapping 
functions in the literature to the genes that they control. 
Unique transcription factor families that regulate genes in 
patterns 2 and 3 include AHRR and EBOX. AHRR, which 
is a unique clock transcription factor family associated 
with pattern 2, has been associated with the regulation of 
cell adhesion and matrix formation [57], which is in agree- 
ment with the functions of genes in pattern 2. EBOX, 
which is associated with the regulation of genes in both 
patterns, has implicit functions in innate immunity, and it 
has been shown that the mechanism by which TNF-a, a 



well known mediator of the inflammatory response, sup- 
presses clock genes is through interference with EBOX 
mediated transcription [58]. EBOX transcription factors 
have also been shown to be key regulators of the circadian 
rhythm, including well known circadian genes and clock 
controlled transcription factors [59]. Thus, both EBOX 
and AHRR not only have functions that are similar to 
those of their respective gene sets, but also have prior as- 
sociation with circadian regulation, which suggests that 
many of the regulators which control innate immunity 
and metabolism within the liver are at least partially under 
the regulation of the circadian rhythm. 

Alternative factors relevant to the regulation of circadian 
genes 

Our promoter analysis of transcription factors is based 
on the assumption that clock TFs are mainly responsible 
for the transcriptional dynamics that are observed over 
24 hours within the liver. However, other factors, includ- 
ing SCN stimuli and food intake in the form of releasing 
glucocorticoids from the adrenal gland, light stimuli, and 
heat, have been known to impact circadian rhythms 
[60]. This is particularly important in the context of 
genes in pattern 4, whose functions are closely related to 
DNA repair and maintenance and which peaks during 
the night, thereby repairing any UV damage that occurs. 

Among transcription factor families identified from 
Table 1, GREF represents a known transcription factor 
responsive element that is stimulated by glucocorti- 
coids [61]. Glucocorticoids are well known molecules 
that are coordinated by the SCN in order to entrain 
peripheral tissues in a circadian manner [62]. Since 
genes in pattern 1 are controlled by this transcription 
factor family, it suggests that many liver metabolic 
functions may also be regulated by the SCN. This fits 
the current paradigm of SCN regulation, since the 
SCN has been associated with the regulation of feeding 
patterns [63]. Therefore, by controlling metabolic en- 
zymes within the liver to be most active at ZT0, the 
SCN can coordinate metabolic functions to coincide 
with feeding patterns. 

Another alternative regulatory mechanism for circa- 
dian rhythms is through light. Although there is no dir- 
ect light input to liver, light may trigger a signaling 
cascade which eventually activates some functional TFs 
in liver. This activation requires the aid of protein ki- 
nases (PKA) [64], which activate CREB and allow it to 
translocate to the nucleus for transcriptional regulation 
[65]. This two-step activation controls all four circadian 
pattern motifs through CREB factors, and may, in 
addition to CLOCK:BMALl, be responsible for the diur- 
nal activity observed within the liver. The importance of 
light in the sleep/wake cycle has been well established 
[66], and thus it is not surprising that the circadian 
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activity of the liver may be entrained to the same light/ 
dark rhythms that govern sleeping cycles in mammals. 

Finally, genes in pattern 2 & 3 are also regulated by 
the HEAT transcription factor family, which is known to 
be responsive to temperature cues [67]. These cues have 
been shown to regulate the circadian rhythm, and thus, 
it is unsurprising that the innate immune response and 
cell repair mechanisms are tied to heat sensitivity. The 
link between heat shock transcription factors on a circa- 
dian level and innate immune response genes may be a 
part of why burn injuries often illicit severe immune re- 
sponses in the rat liver [68]. Overall, these three alterna- 
tive stimuli likely act in parallel with clock transcription 
factors, in order to give the circadian rhythm plasticity 
in its response to various environmental challenges. 
Figure 4 reinforces the fact that environmental factors 
do affect the circadian transcriptional regulation and 
provides a link to reveal how they affect as well as 
which genes are potential targets of which environ- 
mental factors. Further work should focus on the de- 
grees of control that each of these stimuli has on the 
overall response of the circadian pattern, using a sys- 
tems biology approach to elucidate and refine circa- 
dian networks and the ways in which they crosstalk. 

Conclusions 

Microarrays were used to characterize the gene expression 
within rat livers over 24 hours, and various computational 



tools have been applied to identify and distribute differen- 
tially expressed genes from the microarray data into four 
distinct patterns. Pattern 1 is primarily related to central 
metabolism with the peak at circadian time ZTO. Pattern 
2 is related to cell repair and cell junction formations, and 
peaks at circadian time ZT4. Pattern 3 is primarily related 
to the innate immune response and peaks at circadian 
time ZT8, while Pattern 4 is related to DNA repair and 
peaks at circadian time ZT16. The study also proposed a 
unique regulatory structure for circadian rhythms in the 
liver, which aims to predict the ways in which circadian 
patterns can potentially impact injury. By using techniques 
from promoter analysis coupled with the proposed hy- 
pothesis of clock TF regulation, a set of transcription 
factors and putative regulatory structure have been 
established which show to be consistent in agreement 
with functional and temporal annotations. These tran- 
scription factor dynamics appear to play a strong role 
in circadian regulation, and therefore potentially im- 
pact the liver s response to environmental challenges at 
different circadian times. In conclusion, our study aims 
at (1) unsupervised identification of circadian-relevant 
genes and (2) exploration of potential regulatory mecha- 
nisms relevant to the transcriptional regulation of circa- 
dian genes. While doing that, we eventually proposed a 
promising computational framework that can predict rele- 
vant transcriptional factors and putative regulatory struc- 
ture of circadian genes in a genome-wide manner, 
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establishing a platform towards large-scale comparisons of 
the regulation in homeostatic systems across multiple tis- 
sues and/or conditions. 

Methods 

Animal model 

Male Sprague-Dawley rats (Charles River Labs, 
Wilmington, MA) weighing between 150 and 200 g were 
utilized for this study. The animals were housed in a 
temperature-controlled environment (25°C) with a 12-hour 
light-dark cycle and provided water and standard chow ad 
libitum. Because the animals are nocturnal, this establishes 
a 12 hour LD cycle for the animals, where ZT0-ZT12 rep- 
resents the resting portion of the day, while ZT12-ZT0 rep- 
resents the active portion of the day. All experimental 
procedures were carried out in accordance with National 
Research Council guidelines and approved by the Rutgers 
University Animal Care and Facilities Committee. 

Animals are sacrificed at 9 am, 11 am, 1 pm, 5 pm, 
1 am, and 9 am (~0, 2, 4, 8, 16, and 24 hr zeitgeber 
time). Animals were originally used as part of the study 
by Yang et al. [68] to analyze gene expression following 
burn injury, which dictates the choice of time points. 
Ideally, future works would aim at sampling the latter 
parts of the night and day phases equally to observe 
trends throughout the time course. Liver tissues were 
collected and frozen for microarray analysis (n = 3 per 
time point per group). The tissues were lysed and ho- 
mogenized using Trizol, and the RNAs were further 
purified and treated with DNase using RNeasy columns 
(Qiagen). Then cRNAs prepared from the RNAs of liver 
tissues using protocols provided by Affymetrix were uti- 
lized to hybridize Rat Genome 230 2.0 Array (GeneChip, 
Affymetrix) comprised of 31,099 probe sets. 

Microarray data analysis 

Genome expression data analysis included normalization, 
filtering for differential expression, and clustering. DNA chip 
analyzer (dChip) software [69] was used with invariant-set 
normalization and perfect match (PM) model to generate 
expression values. The analysis of variance (ANOVA) test 
implemented in R [70] was applied to filter differentially 
expressed probesets in the dataset with p-value < 0.001. At 
this step, to control the false discovery rate the method pro- 
posed by Storey [71] was used with a positive false discovery 
rate q-value cutoff of 0.01. Due to a limitation of data, such 
strict thresholds were used and thus several well- 
known circadian genes were not present in the filtered 
gene set (see Additional file 1, sheet 'filtering'). Subse- 
quently, consensus clustering [16] was applied in order 
to identify, in an unsupervised manner, coherent clus- 
ters of co-expressed genes with p-value < 0.001. The 
biological relevance of the intrinsic responses was also 
characterized by evaluating the enrichment of the 



corresponding subsets in circadian rhythm specific 
pathways using KEGG database through ARRAYTRACK 
(p-value < 0.05) [72] and characterizing functionally indi- 
vidual genes. 

Promoter extraction and processing 

Promoters of genes including all transcript-relevant al- 
ternative promoters were extracted from a rich database 
of promoter information. If there is no experimentally 
defined length suggested by Genomatix, a default length, 
500 bp upstream and 100 bp downstream of the tran- 
scription start site which is the most enriched region of 
transcription factor binding sites (TFBSs) according to 
Genomatix, is used [73]. In order to accelerate the 
process of identifying putative transcriptional regulators, 
promoters are pre-processed as in [74]. Specifically, 
Matlnspector [75] is applied to scan for PWM matches 
on those promoter sequences using optimal parameters 
from MatBase which ensures that the minimum number 
of matches found in non-regulatory sequences i.e. the 
false positive matches is minimized [73]. Each promoter 
of a gene is then re-modelled to become a list of TFBSs 
ordered by their local positions on the promoter se- 
quences and represented by corresponding TF names 
along with their binding orientations. The conversion 
supports for fast search the presence of a TFBS or a set 
of TFBSs on promoter sequences [74]. 

Prediction of relevant transcriptional regulators 

In order to predict putative TFs relevant to the tran- 
scriptional regulation of circadian-relevant gene sets, we 
explored the basic underlying assumption of compara- 
tive genomics which states that functional regions evolve 
in a constrained fashion and therefore at a lower rate than 
non-functional regions [41,42]. Consequently, orthologous 
promoters of each gene are further extracted. Each pro- 
moter P is now characterized by a set of orthologous pro- 
moters from the same gene of other vertebrate species, if 
available (e.g. Homo sapiens, Macacamulatta, Pantroglo- 
dytes, Musmusculus, Equuscaballus, Canis lupus famil- 
iaris, and Bos Taurus). To be consistent in the search for 
conserved regions on promoter sequences in order to 
identify putative TFBSs, we eliminate those that do not 
consist of more than two orthologous promoters [76]. 
DiAlign TF [77] with default parameters (core similarity: 
0.75, matrix similarity: optimal threshold for each PWM 
suggested from MatBase) was applied to identify con- 
served regions on promoter P. We next apply Matlnspec- 
tor [75] to scan for all physical TFBSs and only keep those 
that locate on conserved regions. Subsequently, given a 
set of genes, TFBSs that are enriched with a common 
threshold (70% in this study) are identified and considered 
as TFs relevant to the transcriptional regulation of the 
gene set. 
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Identification of clock-controlled transcription factors 

Besides TF families composing of clock TFs identified 
from literature (CLOCKiBMALl, PARF and RORA [78]), 
we hypothesized that transcription factors that contain 
CLOCKiBMALl binding sites on their promoters and also 
differentially expressed in the homeostatic system are 
likely to be clock- controlled TFs. After identification of 
differentially expressed genes from the entire dataset, 
Matlnspector [75] was applied to scan for TFBSs present 
on corresponding promoters of each genes. Genes that 
contain CLOCKiBMALl binding sites and are also re- 
ported as transcription factors in MatBase [73] are consid- 
ered as clock-controlled TFs. 

Identification of common sets of clock transcription 
factors 

It has been noticed that TFs in higher eukaryotes regulate 
gene expression in a combinatorial manner rather than in 
isolation [31,32]. Consequently, to explore the underlying 
mechanisms of clock TFs in the relation to the regulation 
of these circadian rhythms, we hypothesize that significant 
sets of clock TFs present on promoter sets of co- 
expressed circadian-relevant genes are more likely to be 
coordinated regulators in modulating the transcriptional 
process of circadian rhythms. We thus applied the breadth 
first search technique to identify set of clock TFs that are 
commonly present on the promoter set of genes in each 
circadian pattern. Since a gene may have multiple alterna- 
tive promoters, if a set of clock TFs present on any pro- 
moter of the gene it is considered as presence on that 
gene. The procedure first identifies all potential clock 
TFBSs that are commonly present on the corresponding 
promoter set and then search for all possible combina- 
tions of those TFBSs (common level of 70%). Statistical 
significance of each common set of clock TFs is evalu- 
ated to select significant sets for further examination 
(p-value < 0.05, q-value < 0.10) (detailed in Additional file 1). 

Statistical significance of sets of clock TFs 

In order to construct the relations of clock TFs to the 
regulation of circadian rhythms, we estimate the signifi- 
cance of common sets of clock TFs identified from the 
corresponding gene set of each transcription pattern vs. 
the background set. The procedure proposed by Bozek 
et al. [79] was followed to extract a GC-matched back- 
ground gene set (-10,000 genes) in order to compensate 
the relatively high GC-content promoters of circadian 
relevant genes (detailed in Additional file 1). 

Given a gene set G, the significance of set of clock TFs 
A is calculated based on a hyper-geometric p-value de- 
fined as follows: 



p-value (A) 




where B and b is the number of genes and the number 
of hits respectively in the background set; N and n is the 
number of genes and hits in gene set G, respectively. A 
'hit' is the presence of set A on a promoter of a gene. 

Additionally, to restrain the multiple testing problem 
in searching for significant sets of clock TFs, the boot- 
strap method proposed by Storey [71] was utilized to es- 
timate the positive false discovery rates (q-value). A rate 
of 10% was used in this case. 

Additional file 



Additional file 1: This file includes all the information pertaining to 
this paper including gene names, raw expression data, binding 
sites, transcription factors and all relevant statistics. 
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